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Abstract. We consider edge detection as the problem of measuring and localizing changes 
of light intensity in the image. As discussed by Torre and Poggio (1984), edge detection, 
when defined in this way, is an ill-posed problem in the sense of Hadamard. 

Using standard regularization theory, we regularize the problem with a stabilizing 
functional that is a specific form of a Tikhonov stabilizer, following Reinsch (1967) and 
Schoenberg (1964). The regularized solution that arises is then the solution to a variational 
principle. In the case of exact data, one of the standard regularization methods (see Poggio 
and Torre, 1984) leads to cubic spline interpolation before differentiation. We show that In 
the case of regularlyrspaced data this solution corresponds to a convolution fiiler--to be 
applied to the signal before differentiation— which is a cubic spline. In the case of non-exact 
data, which is the most interesting situation, we use another regularization method that leads 
to a different variational principle. We prove (1) that this variational principle loads to a 
convolution filter for ;he problem of one-dimensional edge detection, (2) that the form of 
this filter is very similar to the gaussian filter, and (3) that the regularizing parameter X In 
the variational principle effectively controls the scale of the filter. 

Finally, we outline several issues arising from our solution to the edge detection 
problem: (1) the use of methods from regularizing theories for finding the optimal value 
of the regularizing parameter X; (2) the connection between these methods and the scale- 
space method for edge detection; (3) the relationship between our edge detector and 
oliier detectors, especially the Marr/Hi!dreth edge detector; (4) the extension of our one- 
dimensional solution to two-dimensional edge detection; and (5) the extension of our method 
to deal with differentiation of surface data (though the physical constraint underlying the 
form of the regularizer is not valid in general for depth data); this issue is connected to the 
problem of interpolating and approximating depth data. 
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1. Introduction 

Edge detection does not have a precisely defined goal. The word "edge" itself, which refers 
to physical properties of objects, is somewhat of a misnomer. Several years of experience 
have shown that the ideal goal of detecting and locating physical edges in the surfaces 
being imaged is very difficult and still out of reach (for a review see Brady, 1982). Edge 
detection has come to be defined as the first step in this goal of detecting physical changes 
such as object boundaries— the operation of detecting and locating changes in intensity in 
the image. Other processes which operate on these measurements of intensity changes will 
then group boundaries and label and characterize them in terms of the properties of the 
3-D surfaces. 

Intended in this narrow sense, edge detection— this first step in processing the 
image— is mainly the process that measures, detects and localizes changes of intensity. 
Derivatives must be estimated correctly to label the critical points in the image intensity 
array, characterize their local properties (are they minima or maxima or saddle points?) 
and thus relate them to the underlying physical process (are they shadow edges or depth 
discontinuities?). As a consequence, several different derivatives of the image, possibly at 
different scales, may have to be estimated.^ 

In this sense, Torre and Poggio (1984) considered edge detection as a problem of 
numerical differentiation of images. The problem is not straightforward, and attempts over 
many years have proven its difficulties. Considered as a problem of numerical differentiation, 
edge detection turns out to be an ill-posed problem. As explained by Poggio and Torre 
(1984), mathematically ill-posed problems are problems where the solution either does not 
exist or is not unique or does not depend continuously on the data. 

Numerical differentiation is a (mildly) ill-posed problem because its solution does not 
depend continuously on the data. It is therefore natural to try to solve this problem by using 
regularization techniques developed in recent years for dealing with mathematically ill-posed 
problems. The problem can be regularized by the use of a wide class of filters (Torre and 
Poggio, 1984, section 2.4; see also Duda and Hart, 1973). In the following section we 
consider two specific regularizing operators, that in some sense are very natural. 

2. Regularizing Edge Detection 

To regularize an ill-posed problem and make it well-posed, one has to introduce generic 
constraints on the problem. In this way, one attempts to force the solution to lie in a subspace 
of the solution space, where it is well defined. The basic idea of regularization techniques 
is to restrict the space of acceptable solutions by choosing the function that minimizes an 
appropriate functional. Poggio and Torre consider in particular standard regularization 
theory based on quadratic variational principles. They list three main techniques for 
regularizing the ill-posed problem of finding z from the data y such that Az = y. They 
involve the choice of norms ||-|| (usually quadratic) and of a stabilizing functional \\Pz\\. 
The choice is dictated by mathematical considerations, and, most importantly, by a physical 
analysis of the generic constraints on the problem. Three main methods can then be applied 
(see Bertero, 1982): 

(1) Among z that satisfy ||P2|| < C, where (7 is a constant, find z that minimizes 

U^-y\\, (1) 



' A very similar problem arises in the characterization of surface properties— in particular their 
differential properties— from depth data. 
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Figure 1 Cubic Interpolating Spline Filter (L4 function of Schoenberg) 

(2) Among z that satisfy \\Az - y\\ < C, find z that minimizes 

ll^^ll, (2) 

(3) Find z that minimizes 

\\Az-y\\'+\\\Pz\\\ (3) 

where X is a regularization parameter. 

The first method consists of finding the function z that satisfies the constraint ||F2|| < C and 
best approximates the data. The second method computes the function z that is sufficiently 
close to the data {C depends on the estimated errors and is zero if the data are noiseless) 
and is most "regular". In the third method, the regularization parameter X controls the 
compromise between the degree of regularization of the solution and its closeness to the 
data. 

In the case of edge detection considered as numerical differentiation, we want an 
approximation / to the intensity data y, at sample points i, that is well behaved under 
differentiation. Thus we consider an operator A which samples the function / on the lattice 
such that Af\xi = fUi for t = l,...,N. 

The problem is then to find a suitable norm and a suitable stabilizing functional ||F/||. 
It is natural to chose for P the simplest form of Tikhonov's stabilizing functionals (Tikhonov 
and Arsenin, 1976) with P = -£][ and the usual L2 norm. For A; = 2 this choice corresponds 

to a constraint of smoothness on the approximated intensity profile z, with \\Pf\\ = 0. 
Its physical justification is that the noiseless image has to be smooth in the sense that its 
derivatives must be bounded because the image is band-limited by the optics. Band-limited 
functions have bounded derivatives because /' < fiM, where M = sup F{u), Q is the 
cut-off frequency, and F{u) is the Fourier transform of f{x). Physically, the constraint of 
smoothness allows us to effectively eliminate the noise that creeps in, after, or during the 
sampling and transduction process and makes the operation of differentiation unstable. We 
stress that this is not the only stabilizing functional possible for this problem, although it is 
probably the simplest one. 

2.1. The second regularization method and Interpolating cubic splines 

With this choice of P, the second regularization method is: among / such that /(x.) = y, 
find / that minimize //"' dx. A theorem by Schoenberg (see Greville, 1969) shows that 
the solution to this problem is a cubic spline interpolating between the data points. The 
following result is a reformulation of results due to Schoenberg (1946, 1964): 

Theorem 1: The cubic spline function f interpolating evenly-spaced data points 
and minimizing J f"' dx can be obtained by convolving the data points with a 
cubic spline filter which corresponds to the L4 function of Schoenberg. 

A plot of L4 is given in Figure 1. Note that the size of the filter is fixed with respect to 
the sampling lattice. The filter is at every pixel but the central one, where it is 1. Thus L4 
is an interpolating filter that does not perform any smoothing. 
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Figure 2 Filter R derived by regularization principles. 

2.2. The third regularization method, approximating splines and the edge detection 
filter 

The third regularization method leads to the following problem: Find the / that minimizes 

E(^«-/Mf + x/(/"(a:))'^dx (4) 

t 

where X is the regularization parameter which can be found as described later. This problem 
was considered originally by Reinsch (1967) in the case of numerical differentiation and by 
Schoenberg (1964) for the problem of graduation. Both Schoenberg and Reinsch gave 
the solution in terms of approximating cubic splines. In addition, we prove that for most 
practical purposes, the approximating spline function can be obtained by convolving the 
data point y, with the cubic spline convolution filter R shown in Figure 2 (see also Appendix 
1). We then have the following theorem: 

Theorem 2: The solution to equation (4)— the reqularized solution to the problem 
of numerical differentiation — in the case of inexact data, can be obtained by 
convolving the data with a convolution filter which is (a) a cubic spline, and (b) 
very similar to a gaussian. 

Although this result is especially significant in the context of edge detection where 
the search for an optimal filter and its justification has been a longstanding preoccupation, 
it is somewhat surprising that this result does not seem to have been widely appreciated 
in the numerical analysis literature. The exact assumptions under which Theorem 2 is 
valid are discussed in Appendix 1. First, the data must be given on a regular grid (as is 
the case for an image). Second, the image data must either go to zero at infinity or be 
periodic. Under these conditions, the filtering operation is space-invariant and linear (the 
Euler-Lagrange equations corresponding to the quadratic variational problem are linear). 
Thus the approximating spline can be obtained by a convolution operation. Note that the 
result that the regularizing operator corresponding to a quadratic variational principle is a 
convolution filter— for data on a regular grid and toroidal boundary conditions— is valid, In 
general, beyond the case of numerical differentiation. 

Theorem 3: Quadratic, Tikhonov type regularization principles are equivalent to 
convolving the data with a generalized spline filter, if the data are given on an 
regular lattice and the boundary conditions are appropriate. 

A generally interesting question is the physical correctness of the regularized solution. 
In the case considered in this paper the answer is simple and not very Insightful: a necessary 
and sufficient condition for the regularized solution to be correct is that the true intensity 
distribution is a polynomial of order less than 4 between sampling points. This property can 
be derived directly from equation (3) of Appendix 1.2. 

3. Regularization parameter and comparison with the gaussian filter 

Figure 2 shows the filter R obtained by solving the variational principle equation (1) in 
Appendix 1.1. Its shape and size depends on the regularizing parameter X. Figure 3(a) 
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(a) 



(b) 



Figure 3 Filter R and first derivative for different values of regularizing parameter X. (a) 
X affects size of (the first derivative of) R but (b) does not appear to affect shape of 
R. (Amplitudes are normalized in (a); both amplitudes and widths are scaled linearly and 
independently in (b) for comparison.) 




R filter 
Gaussian 




(b) 




(c) 

Figure 4 Comparison of one-dimensional regularizing filter R with Gaussian: (a) zeroth, 
(b) first, and (c) second derivatives. 



shows the first derivative of the filter for different values of X. The continuous version of 
the filter, derived in Appendix 2, is practically indistinguishable from the discrete filter, as 
shown by numerical comparisons. 

It is rather intuitive that the smoothing parameter X controls the effective size of the 
filter. From our numerical work, it seems that X does not significantly affect the shape of 
the filter, but only its size, as shown in Figure 3(b). Changing X amounts to scaling the size 
of the filter up or down. If X is small, smoothness is unimportant, and the filter will tend to 
be an interpolating filter and therefore be similar Xo a 6 function. On the other hand, with 
a very large X, the main weight Is on smoothness, and the filter will tend to be very large. 
The continuous form of the filter suggests that the role of X is indeed equivalent to the role 
of (T for the gaussian (X ~ a'*, as shown in Appendix 4). 

The regularization filter derived here appears to be quite similar to the Gaussian 
distribution. Graphs of the filter R, its first and second derivatives are shown with those 
of the Gaussian in Figure 4. Marr and Hildreth (1980; Hildreth, 1980) have argued that 
the Gaussian is an optimal smoothing filter for detection due to its localization properties 
in both the spatial and frequency domains. The fact that the Gaussian is quite similar to 
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the Optimal filter derived here using regularization principles provides further mathematical 
justification for the use of a Gaussian like filter in edge detection. 

If the boundary conditions are not periodic or natural, then the derivation of Reinsch 
(see Appendix 1.1) provides the correct Green function for those boundary conditions. 
In numerical experiments v^e have found that the Green function obtained in this way is 
space-invariant for all points but those close to the boundaries. 

4. Discussion 

Several questions and extensions suggest themselves in a natural way. Here we list some 
of them. 

4.1. Finding the Optimal X 

Regularization theory can give the optimal value of X if the errors on the smoothing criteria 
and the error of the approximations are known in advance, if the integral of /" is less 
than E, and the sum of (y, - /(x,))^ is less than e, then X = e/E (see Bertero, 1981; 
Tikhonov, 1963; Tikhonov and Arsenin, 1977). Normally, however, errors on the data or 
on the smoothness conditions are not known in advance. Regularization theory provides 
several methods for finding the optimal smoothing parameter X under this circumstance. 
We want to indicate here two main methods: (1) Tikhonov's method, for convolution type 
problems, as is the case here, and (2) the cross-validation method and the generalized 
cross-validation method (Wahba, 1980). 

We plan to evaluate these methods for finding the optimal X for edge detection and 
to test them on real images. An interesting issue that we are also planning to explore is 
the following: The basic idea of the generalized cross-validation method is to check the 
goodness of approximation for each value of X. In order to do that, one computes the 
approximation by using not all but only some of the data points. Thus, the data points that 
are not used for computing the approximations serve as the control points for the goodness 
of the approximation. If one computes the goodness of the fit at different X, one can then 
choose the optimal X. This idea has obvious connections with the use of fingerprints (Yuille 
and Poggio, 1983) for finding the natural scale of the filter (Witkin, 1983); this point is 
discussed next. 

4.2. Optimal X and natural scale 

The size of the filter with which to perform edge detection has always been an unresolved 
issue in computer vision. Our approach makes it clear that one expects, indeed, an optimal 
size of the filter associated with the optimal value of the smoothing parameter X. In more 
recent years, several scales of filtering have been used, partly as a way around this problem. 
Rosenfeld and Kak (1982), Marr and Poggio (1977), Marr and Hildreth (1980; Hildreth, 1980) 
have used several sizes of filters in order to perform edge detection. 

More recently, Witkin (1983) has suggested the use of scale-space filtering, essentially 
filtering across a continuum of scales, as the method by which to choose the optimal scale. 
Witkin suggested some heuristics for picking the natural scale of filtering. We believe that 
cross-validation-type methods may make more rigorous the idea of selecting an optimal 
filtering scale by selecting the optimal X value. We are planning to use fingerprints— the 
zero-crossings across scales of the Laplacian of the convolved image— to find the optimal 
X. 

For scale-space filtering to be maximally effective, the shape of the filter should be 
a gaussian (Yuille and Poggio, 1983a, 1983b; Witkin, 1983). The effect of changing X is 
essentially equivalent to changing the size of the filter, and furthermore, the underlying filter 
is very similar to a gaussian. Therefore, increasing the value of X is equivalent to filtering 
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(a) (b) 

Figure 5 Difference of two-dimensional Gaussians as approximations to (a) Lapiacian of 
gaussian, (b) Lapiacian of two-dimensional regularizing filter R^. The ratio between the 
scales {a) of the two Gaussians is 7. 

the signal with a larger and larger gaussian as required in scale-space filtering. Because 
of the many nice properties of gaussian convolution, this also suggests an efficient way of 
obtaining approximations at increasing X. 

4.3. Relation with other edge detectors 

Because of the close similarity of our cubic spline filter to a gaussian, the edge detector 
that we derive in this paper is very similar to edge detectors proposed previously. Marr 
and Poggio (1977) proposed the difference of two gaussians as an approximation to the 
second derivative of a gaussian. Marr and Hildreth (1980; Hildreth, 1980) have shown that 
the second derivative of a gaussian is, indeed, very close to the difference of gaussians. 
J. Canny's filter (Canny, 1983) is very close to the derivative of a gaussian, and Haralick's 
cubic polynomial interpolant (Haralick, 1982) is again similar to Canny's filter. 

Our derivation justifies the use of a gaussian or a filter very close to a gaussian as 
the best filter for edge detection. Regularization theory yields derivative-of-gaussian filters 
as the optimal filter in a simpler, more general, and, we believe, more rigorous way, than 
previous derivations. In particular, our result makes clear that the quasi-gaussian filter 
regularizes the ill-posed problem of numerical differentiation. The regularizing constraint 
here is that the norm of the derivatives in the noise-free image is small. 

It is interesting that we derive a filter very similar to Canny's, based on simpler and 
more general principles that are not restricted to the optimal detection of step edges. It is 
also interesting to note that the second derivative of the regularization filter, like the second 
derivative of the Gaussian, can be approximated by a difference of Gaussians (although 
not as well). While the second derivative of the Gaussian is best approximated by a space 
constant ratio (the ratio of scales of the two Gaussians) 7 » 1.6 (Marr and Hildreth 1980; 
Hildreth, 1980), increasing the rtitio to 7 « A results in a function vMch better fits the main 
(excitatory) lobe of the regularizing filter, as shown in Figure 5. 

4.4. Extension to two dimensions 

Our approach can be extended to two dimensions in several ways. The most straightforward 
method involves the use of directional derivatives. First order directional derivatives are 
taken along several directions. Each one of them is one-dimensional and can be performed 
according to our one-dimensional edge detector scheme. Depending on the specific goal, 
one may then choose the direction that gives the maximum value of the derivative. This 
corresponds to the non-maximum supression scheme used by Canny (1983). It is also 
equivalent to taking the second directional derivative along the gradient and looking for its 
zeros (Torre and Poggio, 1984). 
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(a) (b) 

Figure 6 Cross section of (a) two-dimensional regularization filter H-z and (b) its Laplacian. 

A second method requires directly formulating the regularization principle in two 
dimensions. Instead of equation (2), one would then have the problem of minimizing 



Y.iyij - /(^o))' + ^jj(i' ff'^^ ^y 



(r>) 



The main problem is the choice of the operator P. If we consider a Tikhonov stabilizer (see 
Poggio and Torre, 1984), then a choice for /' that is smooth enough to allow the use of 
second derivatives of the regularized image is P = V'-^V. This choice is used in Appendix 
3 to derive the filter for the two-dimensional continuous case. The filter is shown in Figure 
6. The choice of the derivative to be used on the filter is a separate, important issue that 
we do not address in this paper. Torre and Poggio (1984) discuss the properties of several 
two-dimensional differential operators, including the second directional derivative along the 
gradient. 

If P is chosen to be the quadratic variation or the square Laplacian, the resulting 
approximations, known as as thin plate splines (Wahba, 1980; Terzopoulos, 1984a), are 
not smooth enough for finding zeros of second derivatives of a function^ as implied by 
Terzopoulos (1984b), It may also be interesting to explore the filters resulting from a 
non-linear functional P. In the case of a non-linear P, one cannot use, in general, the 
standard results about uniqueness and other properties of the solution that are available for 
the quadratic case of Tikhonov's stabilizers, because the functional is no longer convex. 

Clearly, formulations of this type are also relevant for the problem of surface interpola- 
tion and approximation in the sense of Crimson (1982) and Terzopoulos (1984a). In the 
case of sparse data, which they considered, the variational principle does not lead to a 
convolution filter, although it does lead to a standard Green function. On a regular grid it 
leads to a convolution filter similar to the gaussian. As a practical implication, evenly-spaced 
surface data (for example, laser range data) may be interpolated or approximated effectively 
by gaussian convolution. Hence, tasks which involve differentiating surface data, such 
as computing lines of curvature (Brady, Ponce, Yuille and Asada, 1984), could use the 
simpler convolution method to smooth the data. Since Reinsch's method (see Appendix 1.1) 
can deal with boundary conditions different from periodic ones, the corresponding Green 
function can be used to prevent smoothing across depth discontinuities. 

The results of applying the Laplacian of the two-dimensional regularization filter and 
the Laplacian of a gaussian to an image are shown in Figure 7. As expected, due to the 
similarilty of the two filters, both edge detection operators yield similar results. 



"We are indebted to Demetri Terzopoulos for this remark 
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(a) 







(b) 




(c) 



Figure 7 Comparison of Laplacian of two-dimensional regularization filter and Laplacian of 
Gaussian as Edge Detectors: (a) image /, (b) zero-crossings of V'^/e.j*/, (c) zero-crossings 
of VH:2*I. 
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Appendix 1: The regularization filter for discrete data in one dimension 

Below we present two methods for proving that the third regularization method described 
in the text yields a convolution filter in the case of evenly-spaced discrete data in one 
dimension with appropriate boundary conditions. Each derivation includes a method for 
computing the filter, which has the form of a cubic spline. 

1.1. The convolution version of Relnsch's optimal filter 

Reinsch (1967) considered the general problem of interpolating and smoothing a sequence 
of (not necessarily evenly-spaced) data points with uncertainties. Given a sequence of 
points (y.) with uncertainties {Syi), representing values of a function at points (x,), where 
i = 1, 2, . . ., n and xi < X2 < • • • < x„, the problem of finding a smooth function which passes 
near each point (x,, y.) can be formulated as finding the function /(x) that minimizes the 
functional (to be compared with equation 4 in the text) 

Using calculus of variations, a cubic spline function is shown to be the optimal function 
satisfying (1), 

/(x) = a,- + 6,(x - x.) -f c,(x - x.)^ -f di{x - x,•)^ x,- < x < x,+ i, (2) 

for t = 1,2, ...,n- 1, and formulas are presented for calculating from the data (x.) and (y,) 
the coefficients (a.), (6,}, (c,), and (d,). 

We specialize Reinsch's derivation as follows: First, the uncertainty associated with 
each data point is assumed to be the same so that all 6yi = l. Equation (1) above thus 
reduces to equation (4) in the text. Second, the data are assumed to be uniformly spaced, 
so that x,+] - Xi == h for i = 1,2, ...,n- 1. Finally, unlike Reinsch's derivation, we assume 
periodic boundary conditions, i.e., that the solution / is a periodic function over », with 
/(x) = /(x ± kn] for integer k. 

With these specializations we show that each set of coefficients can be calculated by 
multipying the data points y by a constant coefficient matrix: 

a = Ay, b = By, c = Cy, d = Dy, 

where A, B, C, and D are circulant matrices representing the filter and its derivatives. 
Hence, the data can be optimally smoothed (and differentiated) by a convolution operation. 

Using calculus of variations, one obtains from (1) these conditions for the optimal 
function /(x): 

/(x.)+ - /(x,)- = (3a) 

/'(x,)+ - /'(x.)- = (36) 

/"(x,)+ - /"(x,)- = (3c) 

/'"(x,)+ - r(x,)- = -HfM - Vil (3d) 

(using Syi = 1) and 

f""{x) = 0, X, < X < x.+i, (4) 
for t = 1, . . ., n (with xi = x„+i), where 

/(*)(^,)±=lim/(*)(x,±£). 

e— ♦O 
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Equations (3-4) show that f{x) is a cubic spline, having the form of (2); hence, 

/(x,)+ = ai (5a) 

f'{xi)+ = bi (56) 

/"(x.)+ = 2ci (5c) 

/"'(x.)+ = Bdi (5rf) 

and, since all x,- - x,_i = h, 

/(x,)_ = a,_i + bi-ih + Ci-ih^ + d,_i^' (6a) 

/'(xi)_ = 6,_i + 2ci-ih + 3d,_i/i* (66) 

/"(x.)_ = 2c,_i + 6rf,_i/i (6c) 

/'"(x,)_ = 6di_i. (6d) 

Now substitution into (3) yields 

a, - a,_i - 6,_i/i - Ci^ih^ - di_ih^ = (7a) 

6,- - 6,_i - 2ci-ih - Sdi^ih^ = (76) 

2c, - 2c,_i - 6d,_i^ = (7c) 

%di - 6rf._i = -i{f{xi) - Vi). {7d) 

Equations (7) can be manipulated to yield 

bih = (a,+i - a,) - c,/i^ - dife' (8o) 

a,+i - 2a,- + a,_i = ^(c,_i + 4c,- + Ci+i)h^ (86) 

(c,-_i - 2c, + Ci+i)/i2 -^ ^(yi - ai)h^ (Sd) 

Using the notation that: 

/ = the n X n identity matrix [ij^k] where 

J- /l. ifj = fc, i = l,...,n 
^' \o, otherwise 

N = {he nxn "next" matrix [uj^k] where 



p = the n X n "previous" matrix = iV^, 



if y = A; — 1 mod n, j" = 1, 

Otherwise 



we define the matrices 



For example, if n = 4, 

1 0\ 

10 

10 

^0 1^ 



T 

Q 



P-2I + N. 



1 = 



P = 



'4/3 1/3 l/3> 

r_|l/3 4/3 1/3 

' 1/3 4/3 1/3 

a/3 1/3 4/3> 



'0 1 0\ 

10 

1 

.1 0^ 



Q = 



N 



^0 1\ 

10 

10 

^0 1 o/ 




10 
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Defining the vectors a = (a.f , 6 = (6.7if , c = {ah^f, d = [dih'^Y , and y = [y^f, for 
i — 1, . . ., n, equations (8) can be expressed as 

b = {N-I)cL-Q-d (9a) 

Qa = Tc (96) 

d=l{N-I)!, (9c) 

^£=?^(|/-4 (9d) 

These can be simplified to 

& = Ai, b = Bi, iL = C^, d = D^, (10) 
where 

C = S:iQ' + &Tr'Q (lla) 

^ = /-|^QC7 (116) 

D=l{N-I)C (lie) 

B = (iV-/)A-(7-D. (lid) 

Since /, N, and P are circulant matrices, and since the set of circulant matrices is 
closed under matrix addition, multiplication, and inversion, the resulting coefficient matrices 
A, B, C, and D are also circulant, representing filters R{x), R'{x), R"{x), and R"'{x). The fact 
that the filters are derivatives of each other follow from (5); hence, since differentiation is a 
linear operator, the filter R", for example, can be used to both optimally smooth and twice 
differentiate the data. 

1 .2 Another derivation of the regularization filter for discrete data in one dimension 

In this derivation we show that the third regularization method, described by equation (4) 
of the text, yields a convolution filter which is a cubic spline. Standard results from the 
calculus of variations guarantee that our solution has continuous second derivatives. 

Again, the problem is to minimize 



\J[r{x)fdx^Y.i^[x,)-y,)\ 



(1) 



We find the minimum by sending /(x) ^ f{x) + 8f[x) and setting the first variation of (1) to 
zero, 



I r{x) Sf{x)dx + ^ifixi) - yi) 6f{xi) = 0. 



(2) 



This yields the Euler-Lagrange equation 

X/""(x) + fix) J^ 6{x - x.) = 5] y,6{x - x,). (3) 

So far, we have deliberately not specified boundary conditions. For infinite, or toroidal, 
boundary conditions, the function /(x) can be determined in terms of the (y.) by a convolution 
filter if and only if the data points (x.) are evenly spaced. This is because the system is 
then translation-invariant^. We will show this explicitly and give a method for constructing 
the convolution filter. 

The function in (1) is convex, and hence has a unique minimum, so there is a unique 
solution for /(x) in (3). Thus, we only need to see whether a convolution can solve (3). We 



'Boundary conditions other than infinite or toroidal will destroy translation invariance. 
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try a solution 



f{x) = R{x)*y{x) = lR{x-M^)d^ 



(4) 

where "*" denotes convolution, R{x) is a filter, and y{x) = ZiVA^ - x.). We substitute (4) 
into (3) and obtain 

X x; yiR""{x - xi) + Y1Y1 yj-^i'^i - ^M^ - ^.o - E y^^i"^ - ^«) = o r^. 

We compare coefficients of y, and obtain 

\R""{x - xi) + ^ R{xj - Xi)6{x - xj) - 6{x - a:.) = (6) 

3 

If the {xi) are not evenly-spaced, then these equations are inconsistent, and no convolution 
filter exists. If the (x.) are evenly spaced, then the set of equations in (6) reduces to a single 
equation 

\R""{x) + X; R{xj)S{x - xj) - S{x) = (7) 

3 

The solutions to (7) correspond to cubic splines "stitched" together at the points {z,}. 
Let Ri{x] denote the solution in the range x, < x < x,+i. We write 

Ri{x) = a.x^ + fiix^ + 7,x + Si (8) 

The splines are stitched together so that R{x), R'{x), and R"{x) are continuous at the points 
{x,}. From (7), we see that R'" has a discontinuity of -^/?(x.) at x.-. It is straightforward to 
find the relations between the /e.(x) and Ri-i{x) in terms of the parameters a.,/3,-, 7.- and 64. 
This gives 

ttn+i = an- —R{nh) 

nh 
Pn+l = /3n - 2X ^("'') 

Tn+l = 7n - ~^R{nh) 

Sn+l =Sn+ ^-rr^R{nh) 
Da 

where h is the spacing between the lattice points, h — x.+i - x,-, and 

R{nh) = an{nhf + /?n(n/i)2 + 7„(n/i) + *„. (10) 
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Appendix 2: The regularization filter for continuous data in one dimension 

We want to find the function f{x) which minimizes the quadratic functional 

j[f{^)-y[^)fdx + \j[f"[x)fdx. (1) 

where y[x) are the data, given on the continuum. We will show that /(x) can be expressed 
as a convolution of the data y{x] with a suitable filter R{x). 

We obtain the Euler-Lagrange equations for (1) by setting the first variation to zero as 
f{x) ^ f{x) + 6f{x). This gives 

I Sf{x){f{x) - y{x))dx + X y Sf{x)r{x)dx = (2) 

for all variations 6f{x) and hence we have 

X/""(x) + /(x) = y(x). (3) 

The solution to the linear equation (3) is given by 

f{x) = Jr{x,x')y{x')dx (4) 

where r{x,x') is the Green function of (3) obeying 

d* 

^■^<''' ^') + K^. ^0 = H^ - A (5) 

Now (5) is a translation invariant equation and our boundary conditions are periodic or at 
infinity, so r{x, x') is a function of {x - x') only. We write 

r{x, x') = r{x - x'). (6) 

For each different value of X, we have a different equation (5) and hence a different 
Green's function which we denote by R{x,\). Note from (5) that 

R{x,0) = 6{x). (7) 

So in the limit as the confidence in the data is complete, the filter becomes the delta 
function. (It is easy to see all this by taking the Fourier transform of equation 3.) 

It is straightforward to find the Green's function of (5) which vanishes at plus and 
minus infinity. This is given by 

i2(x, X) = -i_e-|x|/v/2X^/^ ^^3 (J±_ _ A (8^ 

2X1/4 Vv^Xi/'t 4; ^""^ 

In Fourier space, the transform of R{x, X) can be obtained directly from (5): 

so we can write 

At X = the Green function goes to a delta function, as we saw in (7). 
Define n by 

X/i* = 1. (11) 

Then, for x > 0, we have 
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and 

Thus the extrema of R{x, n) occur at 

x= — nn, n=l,2. .. (14) 

and at these extrema, R takes the value 

So if n is small (X large), the extrema are at large x and correspond to small values of R. If 
fi is large (X small), the nearest extrema occur at small x and correspond to large values of 
R. The function changes sign many times. 
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Appendix 3: The regularization filter for continuous data in two dimensions 

We can extend the results of the previous section to two dimensions. Our generalization of 
the smoothing function f{f"{x))'~dx is 



// 



V'Vf{x)V'Vf{x)dx, (1) 

Courant and Hilbert (1953) show that the Euler-Lagrange equation is 

V'-^V^VV-O. (2) 

We write the regularization of the two-dimensional smoothing problem in the form 

/ jifU) - y{'^)f<ix -f X I liV'Vf{x)V'Vf{x))dx. (3) 

The Euler-Lagrange equations of the combined system are 

XV^V^VV{x)4/(x)^y(x). (^) 

This equation is translation invariant and so, for boundary conditions at infinity or periodic, 
the solution can be written as a convolution 

f{x) = n^ix) * y{x) (5) 

where 

xv2v2v2/e2(x) + «2U) = %) (6) 

where 6{x) is the Dirac delta function. Again, observe that for small X the filter I^x) tends 
to the delta function. 



To solve (6) we take its Fourier Transform, and find 

1 

XcjG + 1 



TR,{u) = -1--. (7) 



This gives a solution 

We express x and w in polar coordinates 

= [r,(l>), u^-={w,0). (9) 



X 

So 

«oo /•27r 



1 f°° f e+''^co>^(9-4,) 
27r Jq Jo \w^ + 1 



(10) 



Now 

.27r 






I e'^''^'^o-~'^)dO~^irJ„{wr) (Ji) 

where Jo is the zero order Bessel function. This gives 

that is, 71*2 is the Hankel transform of l/(Xu/' + l). This integral can be solved numerically. 
The square Laplacian stabilizer leads to a similar formula, i.e. to 

The corresponding filter is not smooth enough to allow second derivatives to be taken, but 
is sufficient for first derivatives. The integral is found in standard texfs (for example, 
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Gradshteyn and Ryzhik, 1980), yielding 
where k(i{x) is a Thompson function given by 

kci{x) = [Inl - C)hci{x) - ^bcr{x) + f]{-\)' "—^ 't 1 (15) 

k=o 2''=+'-^[(2A;-f l)!]^;^, ^'^ 

where 






where a and /3 are constants. 



:««) 



and C = 0.5772... is Euler's constant. The asymptotic behavior is 

kci{x) ^ y ^e-«* sin(-/?x) (17) 
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Appendix 4: The regularization filter as an approximation to the Gaussian 

In this appendix we show that the regularization filter approximates the gaussian, both in 
one and two dimensions. 

4.1 Comparison of one-dimensional filters 

We show that the one-dimensional regularization filter is an approximate solution to the 
diffusion equation and therefore approximates a gaussian. 

The one-dimensional regularization filter in is given by 



We expand this in a Taylor series in ^ to get 

r iiT. \ 

(2) 



2n/2 



,^(g).o,H 











d'^fi 


-V 










dx^ ~ 


2/2 


and 








ail 1 

^^ 2/2 


+ 0{fi 


which 


satisfy, 


to order 


(/^a:)^ 


the equation 

dx^ 


.ah 



This expansion is valid when nx is small, i.e., when x is small compared to X'/^ In this case 
the first two terms which we denote by /^(x,/i) are a good approximation to the function. 
We calculate 

(3) 



{^) 



(5) 
Thus this function obeys the diffusion equation, 

a'^Ji ail 

'a^-m' ^^^ 

with parameter t == i/^-^ ^ ix*/^ in the region where fix is small. As fi decreases, this 
region gets larger and the region in which the function approximates a Gaussian increases. 

This theoretical analysis supports the numerical results (for discrete data) which show 
that R can be approximated by a Gaussian. Furthermore, recalling that the standard 
deviation a of the Gaussian is given by cr = /2l, the analysis shows that the standard 
deviation of the corresponding Gaussian is X^/^. 

A comparison of R with the gaussian G = (e-==V2<»') can also be done directly in the 
Fourier domain, where TR --= 1/(1 + Xw'*) and TG = e""''*'', as shown in Figure 8. 

4.2 Comparison of two-dimensional filters 

We now consider the two-dimensional case. The regularized filter can be written in terms 
of a Fourier integral 



1 f e-«w£ 
2n J 1 + \u^ 



(7) 
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(a) 



(b) 




Figure 8 Comparison of one-dimensional regularizing filter R with gaussian In Fourier 
domain: (a) zeroth, (b) first, and (c) second derivatives. 

We perform tfie transformation i^i »-* X^/^w to obtain, with /x = :^, 
We expand the exponential in a power series 

/£ f 1 

Thus we have 

Note that the linear term drops out due to asymmetry of the integrand. Keeping the first 
two terms on the right hand side of (10), and denoting this approximation by ih[x,^, we 
calculate 

and 



dR 



~2^J rr^'^"^- 



(12) 



Thus as before, the approximation H2{x, X) satisfies the diffusion equation with t proportional 
to X*/^ The exact function of proportionality can be calculated from (12). 

Again, a direct comparison of the Gaussian with our regularizing filter Ro is done easily 
in the Fourier plane. Both filters are circularly symmetric and therefore depend only on 
the radial frequency w. A comparison in Fourier space of the two-dimensional gaussian 
TG'2 = e-"^'"' and the regularizing filter TR.2 --= l/{\w'^ -i I) is shown in Figure 9. 
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(a) 



(b) 



Figure 9 Comparison of two-dimensional filter R-i with two-dimensional gaussian G^ in 
Fourier plane (a) filters and (b) their Laplacians. 
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